In the following, a method of automatically checking that some given algorithm is indeed an integration method of a certain order is given.
This is only for the purpose of discussing sympy; much better methods are available for proving that certain methods have various properties.
Additional comment: while you are encouraged to use sympy or similar tools for making quick studies of certain problems, please be aware that, as any software package, sympy and similar tools almost certainly have bugs somewhere, and you should always have a way to verify the results that you get from your code.
In [1]:
from IPython.display import display
from sympy.interactive import init_session
init_session()
import sympy as sp
x = sp.Symbol('x')
h = sp.Symbol('h')
f = sp.Function('f')
print(f(x))
In [2]:
def difft(e):
return e.diff(x)*f(x)
x0 = x
x1 = difft(x0)
x2 = difft(x1)
x3 = difft(x2)
bla = x0 + x1*h + x2*h**2/2 + x3*h**3/6
In [3]:
print(sp.latex(bla))
Taylor expansion of an ODE solution, up to order 3 in $h$:
\begin{equation} x(t+h) \approx \frac{h^{3}}{6} \left(f{\left (x \right )} \frac{d^{2}}{d x^{2}} f{\left (x \right )} + \left(\frac{d}{d x} f{\left (x \right )}\right)^{2}\right) f{\left (x \right )} + \frac{h^{2}}{2} f{\left (x \right )} \frac{d}{d x} f{\left (x \right )} + h f{\left (x \right )} + x \end{equation}Note that this was constructed by applying the operator $f(x) \frac{d}{dx}$ "by hand". Since it's only for the first 3 terms, this is acceptable. However, it would be nice to do it for any order.?
In [4]:
t = sp.Symbol('t')
h = sp.Symbol('h')
x = sp.Function('x')
f = sp.Function('f')
print(sp.latex(sp.series(x(t+h), x = h, x0 = 0, n = 4)))
This sympy.series functionality seems to be of use, even though the $\frac{d}{dt}$ operator is not explicit in our ODE case.
As usual, you can access the "full" documentation with help:
In [5]:
help(sp.Expr.series)
Personally, I would not consider this fully informative, especially because there's a parameter that is not discussed at all (the logx thing).
But it does give the information we need.
In [6]:
t = sp.Symbol('t')
h = sp.Symbol('h')
x = sp.Function('x')
f = sp.Function('f')
# What follows is a function that will construct the Taylor expansion
# for an ODE solution (i.e. a series in the timestep, with "explicit"
# coefficients made from the right hand side and its derivatives).
# This can easily be generalized to more dimensions if need be.
def get_ode_series(
x, f, t, h, n = 3):
# when computing the terms of the series, we don't actually
# care about time dependency, since we already know how to
# apply the chain rule.
y = sp.Symbol('xi')
difflist = [y]
for i in range(1, n+1):
difflist.append(
(difflist[-1].diff(y)*f(y)))
# in the Taylor expansion, the temporary "y" is replaced with x(t)
return sum(difflist[i].subs(y, x(t))*(h**i)/sp.factorial(i)
for i in range(n+1))
get_ode_series(x, f, t, h)
Out[6]:
In [7]:
# Once we can construct the Taylor series up to any order,
# we can start playing with ODE solvers.
# As long as the solver is not too complicated, we can use sympy to check
# that a solver has a certain order.
def Euler_spec(x0, f, dt):
return x0 + f(x0)*dt
def get_solver_error(
solver,
x, f, t, h, n = 3):
y = solver(x(t), f, h)
s1 = get_ode_series(x, f, t, h, n)
# sympy will do all the hard work for us when we call the "series"
# method; since y is an expression, this is the exact function for
# which we called help earlier.
s2 = y.series(x = h, x0 = 0, n = n)
return sp.simplify(s1 - s2)
In [8]:
# here's what the error term looks like for the Euler method:
get_solver_error(Euler_spec, x, f, t, h, n = 2)
Out[8]:
In [9]:
# But we're actually only interested in the order of the error term,
# and sympy has a function that will give us that:
sp.Order(get_solver_error(Euler_spec, x, f, t, h, n = 3), (h,0))
Out[9]:
In [10]:
# The Heun method is a second order method.
# and we can now "check"
def Heun_spec(x0, f, dt):
y = x0 + dt*f(x0)
return x0 + (f(x0) + f(y))*dt/2
sp.Order(get_solver_error(Heun_spec, x, f, t, h, n = 4), (h,0))
# By the way, I have no idea why that factor of 2 is carried over;
# I guess we can simply ignore it for now.
Out[10]:
In [11]:
def Taylor2_spec(x0, f, dt):
return x0 + dt*f(x0) + f(x0)*f(x0).diff(x0)*dt**2/2
sp.Order(get_solver_error(Taylor2_spec, x, f, t, h, n = 4), (h,0))
Out[11]:
In [12]:
# one step of classic Runge Kutta has an order 5 error:
def cRK_spec(x0, f, dt):
k1 = f(x0)
k2 = f(x0 + dt*k1/2)
k3 = f(x0 + dt*k2/2)
k4 = f(x0 + dt*k3)
return x0 + (k1 + 2*(k2 + k3) + k4)*dt/6
sp.Order(get_solver_error(cRK_spec, x, f, t, h, n = 7), (h,0))
Out[12]:
In [13]:
# at https://en.wikipedia.org/wiki/List_of_Runge%E2%80%93Kutta_methods#Explicit_methods
# you can find a whole bunch of parameters for Runge Kutta methods.
# we might as well have a specific way of checking their error terms
def explicitRK(x0, f, dt,
c = [0],
b = [1],
a = [[]]):
k = [f(x0)]
for i in range(1,len(c)):
k.append(f(x0 + dt*sum(k[j]*a[i][j]
for j in range(len(a[i])))))
return x0 + dt*sum(k[i]*b[i] for i in range(len(b)))
def get_eRK_error(
params,
x, f, t, h, n = 3):
y = explicitRK(x(t), f, h, params[0], params[1], params[2])
s1 = get_ode_series(x, f, t, h, n)
s2 = y.series(x = h, x0 = 0, n = n)
return sp.simplify(s1 - s2)
# and here's the error for Euler:
sp.Order(get_eRK_error(([0], [1], [[]]),
x, f, t, h, n = 5),
(h,0))
Out[13]:
In [14]:
def get_order(
c = [],
b = [],
a = [],
n = 4):
t = sp.Symbol('t')
h = sp.Symbol('h')
x = sp.Function('x')
f = sp.Function('f')
bla = get_eRK_error((c, b, a),
x, f, t, h,
n = n)
return sp.Order(bla, (h,0))
# Here's the error for the "Explicit midpoint method"
get_order(
c = [sp.Rational(0, 1),
sp.Rational(1, 2)],
b = [sp.Rational(0, 1),
sp.Rational(1, 1)],
a = [[],
[sp.Rational(1, 2)]])
Out[14]:
In [15]:
# Error for "Kutta's third-order method"
get_order(
c = [sp.Rational(0, 1),
sp.Rational(1, 2),
sp.Rational(1, 1)],
b = [sp.Rational(1, 6),
sp.Rational(2, 3),
sp.Rational(1, 6)],
a = [[],
[sp.Rational(1, 2)],
[sp.Rational(-1, 1), sp.Rational(2, 1)]],
n = 6)
Out[15]:
In [3]:
def centered_differences(neighbours = 1):
npoints = neighbours*2+1
x = sp.Symbol('x')
y = sp.Symbol('y')
f = [sp.Symbol('f_{0}'.format(i)) for i in range(npoints)]
t = sum(f[i]*y**i/sp.factorial(i) for i in range(len(f)))
h = sp.Symbol('h')
F = sp.Function('f')
solution = sp.solve([t.subs(y, -i*h) - F(x-i*h)
for i in range(-neighbours, neighbours+1)],
f)
expansion = {}
for k in solution.keys():
expansion[k] = sp.simplify(solution[k].series(
x = h, x0 = 0, n = npoints))
return solution, expansion
In [4]:
centered_differences(2)
Out[4]:
In [5]:
a, b = centered_differences(1)
In [6]:
type(a)
Out[6]:
In [9]:
def finite_differences(n = 1, m = 1):
npoints = n+m+1
x = sp.Symbol('x')
y = sp.Symbol('y')
f = [sp.Symbol('f_{0}'.format(i)) for i in range(npoints)]
t = sum(f[i]*y**i/sp.factorial(i) for i in range(len(f)))
h = sp.Symbol('h')
F = sp.Function('f')
solution = sp.solve([t.subs(y, i*h) - F(x+i*h)
for i in range(-n, m+1)],
f)
expansion = {}
for k in solution.keys():
expansion[k] = sp.simplify(solution[k].series(
x = h, x0 = 0, n = npoints))
return solution, expansion
In [10]:
finite_differences(n = 0, m = 1)
Out[10]:
In [11]:
finite_differences(n = 0, m = 2)
Out[11]:
In [12]:
finite_differences(n = 1, m = 1)
Out[12]:
In [ ]: